/*! 
 ***************************************************************************
 * \file   convection4.cpp
 * Forced Convection-Diffusion in NonUniform 2D meshes
 *
 * \par Description
 * This code solves the next equation:
 * \f[ 
 * \frac{\partial T}{\partial t} + 
 * u_j \frac{\partial T}{\partial x_j} = \nabla^2 T
 * \f]. 
 * \par
 * This equation is solved in a unit square \f$ x, y \in [0,1] \f$. 
 * The velocity is prescribed and has the form:
 * \f[
 * u = -A \cos(\pi y) \sin(\pi \lambda x / L_x );
 * \f]
 * \f[
 * v = \frac{A \lambda}{L_x} \sin(\pi y) \cos(\pi \lambda x / L_x);
 * \f]
 * The boundary conditions are as shown in the next figure:
 *
 * \image html convection.png "Unit square in 2D" width=10cm
 * \image latex convection.eps "Unit square in 2D" width=10cm
 *
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 * \par
   \verbatim
   % make
   % ./convection4 \endverbatim
 * 
 ***************************************************************************
 * \author    Luis M. de la Cruz  [ Sat May 23 12:06:36 CDT 2009 ]
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/NonUniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/NUCDS_CoDi.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;

// Functor to define the point distribution.
class MyFunction {      
  double lambda; // distribution factor 
public:                                    
    MyFunction(double l) : lambda(l) { }   

    double operator() (double epsi) {      
	return ( (exp(lambda * epsi) - 1.0) / 
		 (exp(lambda       ) - 1.0) ); 
    }
};

// Scalar field in one dimension
typedef TunaArray<double,1>::huge ScalarField1D;

// Scalar field in two dimensions
typedef TunaArray<double,2>::huge ScalarField2D;

template<class sfield, class t_malla>
void initialVelocity(sfield &u, sfield &v, t_malla &malla,
		      double A, double l_x, double l)
{
    const int bi = u.lbound(firstDim) + 1, ei = u.ubound(firstDim) - 1,
	bj = u.lbound(secondDim) + 1, ej = u.ubound(secondDim) - 1;

    double A2 = A * l / l_x;

    ScalarField1D x = malla.getCoord(X), y = malla.getCoord(Y);

    for(int i = bi; i <= ei; ++i)
	for(int j = bj; j <= ej; ++j) {
	  u(i,j) = -A * cos (PI * y(j)) * sin (PI * l * x(i) / l_x);    
	  v(i,j) = A2 * sin (PI * y(j)) * cos (PI * l * x(i) / l_x);
	}
}
  

int main() 
{

    double length_x, length_y;
    int num_nodes_x, num_nodes_y;
    double right_wall = -0.5, left_wall = 0.5;
    double tolerance, dt, error = 1.0, residual;
    double lambda, amplitude;
    double lcoef;
    int max_iter, nmax, frequency, iteration = 0;

// Reading data from the input file
    std::ifstream input_cin ("input4");
    input_cin >> length_x
	      >> length_y
	      >> num_nodes_x 
	      >> num_nodes_y 
	      >> dt
	      >> nmax
	      >> tolerance
	      >> max_iter
	      >> lambda
	      >> amplitude
	      >> lcoef
	      >> frequency;
    input_cin.close();

    StructuredMesh<NonUniform<double, 2> > mesh(length_x, num_nodes_x,
						length_y, num_nodes_y);
    int mid_nod_x = num_nodes_x / 2 + 1;
    int mid_nod_y = num_nodes_y / 2 + 1;
//
// x-axis distribution
//
    mesh.calcPoints(X, MyFunction(lcoef), 0, length_x * 0.5, mid_nod_x);            
    mesh.calcPoints(X, MyFunction(-lcoef), length_x * 0.5, length_x, mid_nod_x);
    mesh.calcDeltas(X);
//
// y-axis distribution
//
    mesh.calcPoints(Y, MyFunction(lcoef), 0, length_y * 0.5, mid_nod_y);       
    mesh.calcPoints(Y, MyFunction(-lcoef), length_y * 0.5, length_y, mid_nod_y);
    mesh.calcDeltas(Y);
    mesh.print();
    mesh.writeToFile("./Data4/mesh.dat");

    double dx = mesh.getDelta(X);
    double dy = mesh.getDelta(Y);

    int num_vols_x = mesh.getExtentVolumes(X);
    int num_vols_y = mesh.getExtentVolumes(Y);

    int stop;
#ifdef WITH_GNUPLOT
    GNUplot plotter;
    plotter("set view 0,0");
    plotter("splot \"./Data4/mesh.dat\" w l");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> stop;
#endif

    ScalarField2D T(mesh.getExtentVolumes());
    ScalarField2D u(mesh.getExtentVolumes());
    ScalarField2D v(mesh.getExtentVolumes());
    ScalarField2D us(num_nodes_x, num_vols_y); // staggered u-velocity 
    ScalarField2D vs(num_vols_x, num_nodes_y); // staggered v-velocity

    ScalarField2D Tn(mesh.getExtentNodes());
    ScalarField2D un(mesh.getExtentNodes());
    ScalarField2D vn(mesh.getExtentNodes());

    Range all = Range::all();
    T(T.lbound(firstDim), all) = left_wall;   
    T(T.ubound(firstDim), all) = right_wall;

    initialVelocity(un, vn, mesh, amplitude, length_x, lambda);
    NumUtils::staggerFromNodesX(us, un);
    NumUtils::staggerFromNodesY(vs, vn);

#ifdef WITH_DX
    InOut::writeToFile_DX(un, vn, 0, "./Data4/velc.", mesh);
#endif

    SparseMatrix< Diagonal<double, 2 > > A(num_nodes_x, num_nodes_y); 
    //    DiagonalMatrix<double, 2 > A(num_nodes_x, num_nodes_y); 
    ScalarField2D              b(num_nodes_x, num_nodes_y);

    
    ScalarEquation<NUCDS_CoDi<double,2> > conv_diff(T, A, b);
    conv_diff.setDeltaTime(dt);   
    conv_diff.setNonUniformMesh(mesh);
    conv_diff.setGamma(1.0);
    conv_diff.setNeumann(TOP_WALL);
    conv_diff.setNeumann(BOTTOM_WALL);
    conv_diff.setDirichlet(LEFT_WALL);
    conv_diff.setDirichlet(RIGHT_WALL);
    conv_diff.setUvelocity(us);
    conv_diff.setVvelocity(vs);
    conv_diff.print();

    ofstream error_fp("./Data4/error"), residual_fp("./Data4/residual"), 
      nusselt_fp("./Data4/nusselt");

#ifdef WITH_GNUPLOT
    plotter("set grid");
    plotter("set contour");
    plotter("set cntrparam level discrete 0.4, 0.3, 0.2, 0.1, 0, -0.1, -0.2, -0.3, -0.4");
    plotter("set size 0.73,1");
    plotter("unset surface");
#endif

    while ( (error > tolerance) && (iteration <= nmax) ) {
        conv_diff.calcCoefficients();
	Solver::TDMA2DY(conv_diff, tolerance, max_iter);
	residual = conv_diff.calcResidual();
	error = conv_diff.calcErrorL2();
	conv_diff.update();
        nusselt_fp << iteration  << "\t"
		   << NumUtils::calcNusseltVert(T,dx*0.5,dy,0) << endl;
	
	error_fp << iteration  << "\t" << error << endl;
	residual_fp << iteration  << "\t" << residual << endl;

	cout << "\n ---- Iteration = " << ++iteration
	     << " Residual = " << residual
	     << " Error = " << error
	     << endl;

	if (!(iteration % frequency)) {
	  NumUtils::interpolateToNodes(Tn, T, mesh);

#ifdef WITH_DX
	  InOut::writeToFile_DX(Tn, iteration, "./Data4/temp.", mesh);
#endif

#ifdef WITH_GNUPLOT
	  InOut::writeToFile(Tn, 1, "./Data4/temp.", dx, dy);
	  plotter("splot \"Data4/temp.00001\" w l");
	  std::cout << "\n >---> Press <1> and then <enter> to continue \n";	 
	  std::cin >> stop;	
#endif

	}	
    }
#ifdef WITH_DX
    std::cout << "\n\n NOTE: The data were written in the [./Data4] directory in DX format";
#endif

    return 0;
}







